#include <iostream>
#include <fstream>
#include <string>

#include "definitions.hpp"

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>

#include <boost/numeric/ublas/operations.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/triangular.hpp>


#include <boost/numeric/ublas/io.hpp>

Matrix A;
Matrix E;
Vector F;
int N;

bool ReadFile( std::string file);

real_t tau( Matrix & A, Vector & R);
Vector   R( Matrix & A, Vector & U, Vector & F);

int main ( int argc, char *argv[] )
{
    using std::cout;
    using std::endl;
    using boost::numeric::ublas::prod;
    using boost::numeric::ublas::inner_prod;
    using boost::numeric::ublas::norm_1;

    ReadFile( "matrix.txt");
    cout << "N=\t" << N << endl;
    cout << "A:\t" << A << endl;
    cout << "F:\t" << F << endl;
////////////////////////////////////

// Нулевое приблежение и вектор невязки

    Vector U0( N);
    for ( int i = 0; i < N; i++)
    {
        U0( i) = 0.0;
    }

    Vector R0( N);
    R0 = R( A, U0, F);

    cout << "0\t" << U0 << endl;
// Первое приблежение
    Vector U1( N);
    real_t tau_1 = tau( A, R0);

    U1 = prod(  E - A * tau_1, U0) + F * tau_1;

   cout << "1\t" << U1 << endl;
// Последущие итерации
    Vector U_k_p1( N);
    Vector U_k_m1( N);
    Vector U_k( N);
    Vector R_k( N);
    Vector R_k_m1( N);

    real_t tau_k_p1, tau_k, alpha_k_p1, alpha_k, norm_r;
    U_k_m1 = U0;
    U_k    = U1;

    R_k_m1    = R( A, U_k_m1, F);
    R_k       = R( A, U_k   , F);

    tau_k   = tau( A, R_k);
    alpha_k = 1;

    norm_r = norm_1( R_k);

    int k = 1;
    while( norm_r > _PRECISION)
    {
        tau_k_p1 = tau( A, R_k);
        alpha_k_p1 = 1.0 / ( 1.0 - 1.0 / alpha_k * tau_k_p1 / tau_k * inner_prod( R_k, R_k) / inner_prod( R_k_m1, R_k_m1));


        U_k_p1 = prod (( E - A * tau_k_p1) * alpha_k_p1, U_k) + U_k_m1 * ( 1.0 - alpha_k_p1) + F * alpha_k_p1 * tau_k_p1;
        ///
        U_k_m1 = U_k;
        U_k    = U_k_p1;

        alpha_k = alpha_k_p1;
        tau_k   = tau_k_p1;

        R_k_m1 = R_k;
        R_k    = R ( A, U_k_p1, F);

        norm_r = norm_1( R_k);
        //
        //
        cout << k + 1 << ":\t" << U_k_p1 << "\t\tNorm R: " << norm_r  <<  endl;
        cout.flush();
        //
        k++;
    }

    return 0;
}



bool ReadFile( std::string file)
{
    std::fstream fin;

    fin.open( file.c_str(), std::fstream::in);

    fin >> N;

    F.resize( N, false);
    A.resize( N, N, false);
    E.resize( N, N, false);

    for ( int i = 0; i < N; i++)
    {
        real_t tmp;
        fin >> tmp;

        F( i) = tmp;
    }

    for ( int i = 0; i < N; i++)
    {
        for ( int j = 0; j < N; j++)
        {
            real_t tmp;
            fin >> tmp;
            A( i, j) = tmp;
            if ( i == j)
            {
                E( i, j) = 1.0;
            }
        }
    }
    fin.close();
    return true;
}

real_t tau( Matrix & A, Vector & R)
{
    using boost::numeric::ublas::inner_prod;
    return ( - ( inner_prod( R, R) / inner_prod( prod( A, R), R)));
}

Vector   R( Matrix & A, Vector & U, Vector & F)
{
    using boost::numeric::ublas::prod;
    return prod( A, U) - F;
}

